function [D X Di Xi] = ResultsSEP(outputOfAnalyze, subjectNr, electrode, period, doSmoothing)
% subject = outputOfAnalyze.subjects(subjectNr);

if ~exist('doSmoothing', 'var')
    doSmoothing = false;
end

blocks = outputOfAnalyze.results(subjectNr).period{period}.block;
X = []; D = [];
for b = 1:length(blocks)
    if isempty(blocks{b}); continue; end;
    load(blocks{b});
    data = zeros(8,length(results.(electrode).trial)); % N1time, N1value, P1... ... ...
    for t = 1:length(results.(electrode).trial)
        data(1,t) = results.(electrode).trial(t).N1.time;
        data(2,t) = results.(electrode).trial(t).N1.value;
        data(3,t) = results.(electrode).trial(t).P1.time;
        data(4,t) = results.(electrode).trial(t).P1.value;
        data(5,t) = results.(electrode).trial(t).N2.time;
        data(6,t) = results.(electrode).trial(t).N2.value;
        data(7,t) = results.(electrode).trial(t).P2.time;
        data(8,t) = results.(electrode).trial(t).P2.value;
    end
    
    % Retrieve peak info out of the grand average peaksDB
    %[n1 n2 p1 p2 n1v n2v p1v p2v] = RetrievePeakInfo([], subjectNr, electrode);
    [N1 N2 P1 P2 N1Val N2Val P1Val P2Val] = RetrievePeakInfo([], outputOfAnalyze.subjects(subjectNr), electrode);
    % matrix with the same dimensions as data,
    % a cell in this matrix will be the same as its counterpart in data, but then the GA peak is subtracted from this value (retrievesPeakInfo gets these GA values from the peaksDB). 
    
    subGA = zeros(8,size(data,2));
    for t=1:size(data,2)
       subGA(1,t) = data(1,t) - N1;
       subGA(2,t) = data(2,t) - N1Val;
       subGA(3,t) = data(3,t) - P1;
       subGA(4,t) = data(4,t) - P1Val;
       subGA(5,t) = data(5,t) - N2;
       subGA(6,t) = data(6,t) - N2Val;
       subGA(7,t) = data(7,t) - P2;
       subGA(8,t) = data(8,t) - P2Val;
    end
    
    [Db Xb] = NormalizeBlock(subGA);  %[Db Xb] = NormalizeBlock(data);
    D = [D Db];
    X = [X Xb];
    

    [DiB(1,:,b) Xi] = Interpolate(Db(1,:), Xb);
    DiB(2,:,b) = Interpolate(Db(2,:), Xb);
    DiB(3,:,b) = Interpolate(Db(3,:), Xb);
    DiB(4,:,b) = Interpolate(Db(4,:), Xb);
    DiB(5,:,b) = Interpolate(Db(5,:), Xb);
    DiB(6,:,b) = Interpolate(Db(6,:), Xb);
    DiB(7,:,b) = Interpolate(Db(7,:), Xb);
    DiB(8,:,b) = Interpolate(Db(8,:), Xb);

end
[sy si] = sort(X);
D = D(:,si);
X = X(si);

Di = mean(DiB, 3);
if doSmoothing
    
    window = 40;
    for d= 1:8;
        % Roel truukje 
        temp = [repmat(mean(Di(d,1:window)),1,window) Di(d,:) repmat(mean(Di(d,end-window:end)),1,window)];
        temp = fastsmooth(temp(:), window, 3, 1);
        Di(d,:) = temp(window+1:end-window);
  
%         Di(a,:) = fastsmooth(Di(a,:), 40, 3, 1);
    end;
end;

if doSmoothing %Deze doet ie nu niet meer geloof ik
    % % %    instead of smoothing over trials this code smoothes over time (a percentage of the amount of trials, so that when periods of different
    % % %    subjects are added together the smoothing will have been done on the
    % % %    same scale (for example some people may only have had 100 trials, and others
    % % %    500, when smoothing with a predetermined window of 30 trials the
    % % %    window will be comparitively larger for the period of the person
    % % %    with less trials than for the period with more trials.
    % % %    this unequal distribution of smoothing is not ideal when adding
    % % %    together periods with different amount of trials for averaging).

    % New code for variable windowlengths (windowlength now depends on
    % the amount of trials).
    smoothingFactor = 0.1; 
    nrOfTrials = size(X,2);
    window = round(nrOfTrials * smoothingFactor);

    %window = 30; % TODO: Use above instead of this...
%     window = (length(X));

%         if window > (size(X,2)-1)
%             window = (size(X,2)-1);
%         end
    disp('    Smoothing results per s/e/p.');
    for d=1:8
        %includes quick & dirty fix to remove artifact at end of line
        temp = [repmat(mean(D(d,1:window)),1,window) D(d,:) repmat(mean(D(d,end-window:end)),1,window)];
        temp = fastsmooth(temp(:), window, 3, 1);
        D(d,:) = temp(window+1:end-window);
    end
end

end